#age
meanage #already computed
[1] 17.62903
sdage #already computed
[1] 5.763747
#gender
widedf %>% group_by(Gender) %>% tally()
| Gender | n |
|---|---|
| F | 32 |
| M | 30 |
#is there a relationship between age and WASI?
WASIage <- lm(WASI~Age_Z, data = widedf)
WASIagesq <- lm(WASI~Age_Z + I(Age_Z^2), data = widedf)
lmReport(WASIage)
| WASI | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 109.71 | 106.51 – 112.91 | 68.68 | 60.00 | <0.001 |
| Age_Z | -2.07 | -5.29 – 1.15 | -1.29 | 60.00 | 0.203 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.027 / 0.011 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.03 | 0.95 | 0 | 0.18 |
lmReport(WASIagesq)
| WASI | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 108.82 | 104.00 – 113.64 | 45.19 | 59.00 | <0.001 |
| Age_Z | -2.02 | -5.27 – 1.23 | -1.25 | 59.00 | 0.218 |
| Age_Z^2 | 0.91 | -2.74 – 4.55 | 0.50 | 59.00 | 0.620 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.031 / -0.002 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.03 | 0.95 | 0 | 0.18 |
| I(Age_Z^2) | 0.00 | 0.95 | 0 | 0.10 |
#no
#find the number of trials each participant completed by subtracting the number of trials each participant missed
numtrials <- MemDF %>% group_by(SubjectNumber) %>% tally()
numtrials$missed <- 183 -numtrials$n
sum(numtrials$missed)
[1] 37
max(numtrials$missed)
[1] 7
#make a variable defining the correct answer on test trials (2 is correct for FullTrialtype = 9, 11, 12, 13, 14, 15; 1 is correct for FullTrialType = 10)
MemDF$TestCorrect <- ifelse(MemDF$TrialType == 3 & !MemDF$FullTrialType == 10, 2,ifelse(MemDF$TrialType == 3 & MemDF$FullTrialType == 10, 1, NA))
#make accuracy variable (does the response equal the correct response?)
MemDF$TestAccuracy <- ifelse(MemDF$UnshuffledResp == MemDF$TestCorrect & !is.na(MemDF$TestCorrect), 1, ifelse(MemDF$UnshuffledResp == 3-MemDF$TestCorrect & !is.na(MemDF$TestCorrect), 0, NA))
#does test trial performance improve across the task and interact with age?
lmAgeTrialintx <- glmer(TestAccuracy~TrialNum*Age_Z+(1+TrialNum|SubjectNumber), family=binomial, data=MemDF,control = glmerControl(optimizer = "bobyqa", optCtrl=list(maxfun=1e6)))
#quadratic age
lmAgesqTrialintx <- glmer(TestAccuracy~TrialNum*I(Age_Z^2) +Age_Z+(1+TrialNum|SubjectNumber), family=binomial, data=MemDF,control = glmerControl(optimizer = "bobyqa"))
#is quadratic better than linear?
anova(lmAgeTrialintx,lmAgesqTrialintx)
| npar | AIC | BIC | logLik | deviance | Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| lmAgeTrialintx | 7 | 2636.574 | 2677.595 | -1311.287 | 2622.574 | NA | NA | NA |
| lmAgesqTrialintx | 8 | 2638.553 | 2685.434 | -1311.276 | 2622.553 | 0.0208786 | 1 | 0.8851102 |
#quadratic not better, reporting linear
glmerReport(lmAgeTrialintx)
| TestAccuracy | ||||
|---|---|---|---|---|
| Predictors | z | p | Odds Ratios | CI |
| (Intercept) | 13.21 | <0.001 | 4.28 | 3.45 – 5.31 |
| TrialNum | 8.56 | <0.001 | 2.03 | 1.72 – 2.38 |
| Age_Z | 0.51 | 0.612 | 1.06 | 0.86 – 1.30 |
| TrialNum * Age_Z | 0.22 | 0.830 | 1.02 | 0.87 – 1.19 |
| Random Effects | ||||
| σ2 | 3.29 | |||
| τ00 SubjectNumber | 0.51 | |||
| τ11 SubjectNumber.TrialNum | 0.20 | |||
| ρ01 SubjectNumber | 0.91 | |||
| ICC | 0.18 | |||
| N SubjectNumber | 62 | |||
| Observations | 2592 | |||
| Marginal R2 / Conditional R2 | 0.113 / 0.269 | |||
#block performance means
#Block1
pander(mean(widedf$MeanTestBlock1))
0.6325
#Block2
pander(mean(widedf$MeanTestBlock2))
0.8041
#Block3
pander(mean(widedf$MeanTestBlock3))
0.8445
#accuracy for responses to explicit learning questions
#Question: Did this machine always give you the same number of points, or did it sometimes give 0 points and sometimes give you more points?
#get the mean
mean(widedf$RiskySafeAcc)
[1] 0.8451613
#did responses vary by age?
agePredRSAcc <- lm(RiskySafeAcc ~ Age_Z, data = widedf)
#quadratic age
agesqPredRSAcc <- lm(RiskySafeAcc ~ Age_Z + I(Age_Z^2), data = widedf)
anova(agePredRSAcc,agesqPredRSAcc)
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 60 | 1.849635 | NA | NA | NA | NA |
| 59 | 1.843620 | 1 | 0.006015 | 0.1924935 | 0.6624519 |
#quadratic not better, reporting linear
lmReport(agePredRSAcc)
| RiskySafeAcc | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 0.85 | 0.80 – 0.89 | 37.90 | 60.00 | <0.001 |
| Age_Z | -0.02 | -0.06 – 0.03 | -0.88 | 60.00 | 0.382 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.013 / -0.004 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.01 | 0.95 | 0 | 0.13 |
#Question: How many points did this machine give you each time you chose it? (if sub responded that it was safe) OR How many points did this machine give you when it did not give 0 points?
#get the mean
mean(widedf$ValueAcc)
[1] 0.8419355
#does this vary by age?
agePredValAcc <- lm(ValueAcc ~ Age_Z, data = widedf)
#quadratic age
agesqPredValAcc <- lm(ValueAcc ~ Age_Z + I(Age_Z^2), data = widedf)
anova(agePredValAcc,agesqPredValAcc)
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 60 | 2.632215 | NA | NA | NA | NA |
| 59 | 2.596494 | 1 | 0.0357208 | 0.8116825 | 0.3712863 |
#quadratic not better
lmReport(agePredValAcc)
| ValueAcc | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 0.84 | 0.79 – 0.90 | 31.65 | 60.00 | <0.001 |
| Age_Z | 0.02 | -0.04 – 0.07 | 0.65 | 60.00 | 0.516 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.007 / -0.009 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.01 | 0.95 | 0 | 0.11 |
#Remove RTs < 200 ms
MemDF$RT200 <- ifelse(MemDF$RT>.2, MemDF$RT, NA)
#how many short RTs were removed?
MemDF$RTshort <- ifelse(MemDF$RT>.2,0, 1)
sum(MemDF$RTshort)
[1] 22
temp <- MemDF[,c("SubjectNumber","RTshort")] %>% group_by(SubjectNumber) %>% summarize(shortTrials = sum(RTshort))
## `summarise()` ungrouping output (override with `.groups` argument)
#what's the max number of RTs removed for a participant?
max(temp$shortTrials)
[1] 9
#log-transform RTs
MemDF$LogRT <- log(MemDF$RT200)
#trial by age interaction
lmerLogRTAgextrial <- lmer(LogRT ~ TrialNum * Age_Z + (1+TrialNum|SubjectNumber), data = MemDF, control = lmerControl(optimizer = "bobyqa", optCtrl=list(maxfun=1e6)))
#trial by age squared
lmerLogRTAgesqxtrial <- lmer(LogRT ~ TrialNum * I(Age_Z^2) +Age_Z + (1+TrialNum|SubjectNumber), data = MemDF, control = lmerControl(optimizer = "bobyqa", optCtrl=list(maxfun=1e6)))
anova(lmerLogRTAgextrial,lmerLogRTAgesqxtrial)
## refitting model(s) with ML (instead of REML)
| npar | AIC | BIC | logLik | deviance | Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| lmerLogRTAgextrial | 8 | 11208.17 | 11266.82 | -5596.086 | 11192.17 | NA | NA | NA |
| lmerLogRTAgesqxtrial | 9 | 11211.11 | 11277.09 | -5596.553 | 11193.11 | 0 | 1 | 1 |
#age squared not better
lmerReport(lmerLogRTAgextrial)
| LogRT | |||||
|---|---|---|---|---|---|
| Predictors | t | df | p | Estimates | CI |
| (Intercept) | -8.44 | 11279.00 | <0.001 | -0.32 | -0.40 – -0.25 |
| TrialNum | 1.14 | 11279.00 | 0.254 | 0.01 | -0.01 – 0.03 |
| Age_Z | -1.07 | 11279.00 | 0.283 | -0.04 | -0.12 – 0.03 |
| TrialNum * Age_Z | -2.10 | 11279.00 | 0.036 | -0.02 | -0.04 – -0.00 |
| Random Effects | |||||
| σ2 | 0.15 | ||||
| τ00 SubjectNumber | 0.09 | ||||
| τ11 SubjectNumber.TrialNum | 0.00 | ||||
| ρ01 SubjectNumber | -0.17 | ||||
| ICC | 0.38 | ||||
| N SubjectNumber | 62 | ||||
| Observations | 11287 | ||||
| Marginal R2 / Conditional R2 | 0.009 / 0.388 | ||||
#equal expected value trials
mean(widedf$MeanRiskOverallEqEV)
[1] 0.3728879
sd(widedf$MeanRiskOverallEqEV)
[1] 0.2054037
#was risk taking higher or lower than .5
riskt <- t.test(widedf$MeanRiskOverallEqEV, mu = .5)
pander(riskt)
| Test statistic | df | P value | Alternative hypothesis | mean of x |
|---|---|---|---|---|
| -4.873 | 61 | 8.181e-06 * * * | two.sided | 0.3729 |
rstatix::cohens_d(data = widedf, MeanRiskOverallEqEV ~ 1, mu = .5, ci = TRUE, conf.level = .95)
| .y. | group1 | group2 | effsize | n | conf.low | conf.high | magnitude |
|---|---|---|---|---|---|---|---|
| MeanRiskOverallEqEV | 1 | null model | -0.6188405 | 62 | -0.95 | -0.37 | moderate |
#linear age
lmriskageEq <- lm(MeanRiskOverallEqEV~Age_Z, data = widedf)
#quadratic age
lmriskagesqEq <- lm(MeanRiskOverallEqEV~Age_Z+I(Age_Z^2), data = widedf)
anova(lmriskageEq,lmriskagesqEq)
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 60 | 2.560363 | NA | NA | NA | NA |
| 59 | 2.375800 | 1 | 0.1845629 | 4.583385 | 0.0364259 |
#significantly better
lmReport(lmriskagesqEq)
| MeanRiskOverallEqEV | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 0.31 | 0.24 – 0.39 | 8.17 | 59.00 | <0.001 |
| Age_Z | -0.01 | -0.06 – 0.04 | -0.44 | 59.00 | 0.662 |
| Age_Z^2 | 0.06 | 0.00 – 0.12 | 2.14 | 59.00 | 0.036 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.077 / 0.046 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.01 | 0.95 | 0 | 0.11 |
| I(Age_Z^2) | 0.08 | 0.95 | 0 | 0.29 |
lmriskByAgeEq <- twolines(MeanRiskOverallEqEV ~ Age,data=widedf)
ggplot(data=widedf, aes(x=Age, y=MeanRiskOverallEqEV)) +
geom_point() +
stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black") +
scale_y_continuous(breaks=c(0,.25,.5,.75,1), labels = c("0%","25%","50%","75%","100%"),limits = c(0, 1)) +
ylab("Percent Probabilistic Choices \nEqual-EV Trials") +
geom_hline(yintercept=0.5, linetype="dashed") +
theme(text=element_text(family="sans",size=16))
#linear age
lmriskageUnEq <- lm(MeanRiskOverallUnEqEV~Age_Z, data = widedf)
#quadratic age
lmriskagesqUnEq <- lm(MeanRiskOverallUnEqEV~Age_Z+I(Age_Z^2), data = widedf)
anova(lmriskageUnEq,lmriskagesqUnEq)
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 60 | 4.603064 | NA | NA | NA | NA |
| 59 | 4.252430 | 1 | 0.3506334 | 4.864835 | 0.0313176 |
#significantly better
lmReport(lmriskagesqUnEq)
| MeanRiskOverallUnEqEV | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 0.46 | 0.36 – 0.56 | 9.06 | 59.00 | <0.001 |
| Age_Z | -0.04 | -0.10 – 0.03 | -1.03 | 59.00 | 0.307 |
| Age_Z^2 | 0.09 | 0.01 – 0.16 | 2.21 | 59.00 | 0.031 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.096 / 0.065 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.02 | 0.95 | 0 | 0.17 |
| I(Age_Z^2) | 0.08 | 0.95 | 0 | 0.30 |
lmriskByAgeUnEq <- twolines(MeanRiskOverallUnEqEV ~ Age,data=widedf)
ggplot(data=widedf, aes(x=Age, y=MeanRiskOverallUnEqEV)) + geom_point() +
stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black") +
ylab("Percent Probabilistic Choices \nUnequal-EV Trials") +
geom_hline(yintercept=0.5, linetype="dashed") +
theme(text=element_text(family="sans",size=16))
#linear age
lmriskage <- lm(MeanRiskAll~Age_Z, data = widedf)
#quadratic age
lmriskagesq <- lm(MeanRiskAll~Age_Z+I(Age_Z^2), data = widedf)
anova(lmriskage,lmriskagesq)
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 60 | 2.728566 | NA | NA | NA | NA |
| 59 | 2.489727 | 1 | 0.2388389 | 5.659856 | 0.0206117 |
#significantly better
lmReport(lmriskagesq)
| MeanRiskAll | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 0.37 | 0.29 – 0.44 | 9.38 | 59.00 | <0.001 |
| Age_Z | -0.02 | -0.07 – 0.03 | -0.76 | 59.00 | 0.449 |
| Age_Z^2 | 0.07 | 0.01 – 0.13 | 2.38 | 59.00 | 0.021 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.099 / 0.069 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.01 | 0.95 | 0 | 0.14 |
| I(Age_Z^2) | 0.10 | 0.95 | 0 | 0.32 |
lmriskByAge <- twolines(MeanRiskAll ~ Age,data=widedf)
ggplot(data=widedf, aes(x=Age, y=MeanRiskAll)) + geom_point() +
stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black") +
stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black") +
ylab("Percent Probabilistic Choices Overall") +
geom_hline(yintercept=0.5, linetype="dashed") +
theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank())
#median model AIC
median(widedf$TD_AIC)
[1] 139.9945
median(widedf$TDRS_AIC)
[1] 123.8801
#plot TD model AIC by age
ggplot(data = widedf, aes(x = Age, y = TD_AIC)) + geom_point()
#Age regression
TDmodFitAge <- lm(TD_AIC~Age, data = widedf)
lmReport(TDmodFitAge)
| TD_AIC | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 132.91 | 119.53 – 146.28 | 19.88 | 60.00 | <0.001 |
| Age | 0.15 | -0.57 – 0.87 | 0.41 | 60.00 | 0.681 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.003 / -0.014 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age | 0 | 0.95 | 0 | 0.09 |
TDmodFitAgeSq <- lm(TD_AIC~Age + I(Age^2), data = widedf)
lmReport(TDmodFitAgeSq)
| TD_AIC | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 101.42 | 61.39 – 141.46 | 5.07 | 59.00 | <0.001 |
| Age | 4.19 | -0.71 – 9.09 | 1.71 | 59.00 | 0.093 |
| Age^2 | -0.12 | -0.25 – 0.02 | -1.67 | 59.00 | 0.101 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.048 / 0.015 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age | 0.00 | 0.95 | 0 | 0.09 |
| I(Age^2) | 0.05 | 0.95 | 0 | 0.22 |
#no significant age pattern
#Plot TDRS model AIC by age
ggplot(data = widedf, aes(x = Age, y = TDRS_AIC)) + geom_point()
#Age Regression
TDRSmodFitAge <- lm(TDRS_AIC~Age, data = widedf)
lmReport(TDRSmodFitAge)
| TDRS_AIC | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 128.57 | 108.92 – 148.22 | 13.09 | 60.00 | <0.001 |
| Age | -0.50 | -1.56 – 0.56 | -0.94 | 60.00 | 0.349 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.015 / -0.002 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age | 0.01 | 0.95 | 0 | 0.14 |
TDRSmodFitAgeSq <- lm(TDRS_AIC~Age + I(Age^2), data = widedf)
lmReport(TDRSmodFitAgeSq)
| TDRS_AIC | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 146.82 | 86.83 – 206.82 | 4.90 | 59.00 | <0.001 |
| Age | -2.84 | -10.18 – 4.50 | -0.77 | 59.00 | 0.442 |
| Age^2 | 0.07 | -0.14 – 0.27 | 0.64 | 59.00 | 0.521 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.022 / -0.012 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age | 0.01 | 0.95 | 0 | 0.14 |
| I(Age^2) | 0.01 | 0.95 | 0 | 0.11 |
#no significant age pattern
#make a boxplot showing the relationship between AIC for the two models
ggpaired(widedf,"TD_AIC","TDRS_AIC",line.color = "#7B1979", line.size = 0.4,ggtheme=theme_bw()) +
ylab("Model Fit (AIC)") +
xlab("Model") +
scale_x_discrete(labels = c('TD','RSTD')) +
theme(text=element_text(family="sans",size=16),legend.position="none",panel.grid.major = element_blank(), panel.grid.minor = element_blank())
#how many people ar better fit by each model?
widedf$TDRS_Better <- factor(ifelse(widedf$TDRS_AIC<widedf$TD_AIC, 1, 0))
#TDRS Better
sum(as.numeric(widedf$TDRS_Better) == 2)
[1] 42
#TD Better
sum(as.numeric(widedf$TDRS_Better) == 1)
[1] 20
#code for AI magnitude and valence
widedf$AsymmIdxMag <- abs(widedf$TDRS_AsymmIdx)
widedf$AsymmIdxVal <- factor(ifelse(widedf$TDRS_AsymmIdx>0, 1,0))
levels(widedf$AsymmIdxVal) <- c("Positive","Negative")
#plot which model is a better fit as a function of AI valence and magnitude
ggplot(data = widedf, aes(x = TDRS_Better, y = AsymmIdxMag)) +
geom_jitter(aes(color = AsymmIdxVal)) +
theme_bw() +
scale_x_discrete(labels = c('1-learning-rate\nmodel fits better','2-learning-rate\nmodel fits better')) +
scale_color_manual(values = c("#7B1979","#797979")) +
theme(text=element_text(family="sans"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
labs(x = "Model", y= "Asymmetry Index\nMagnitude",
color = "Asymmetry Index\nValence")
widedf$TDRS_TD_AIC <- widedf$TDRS_AIC-widedf$TD_AIC
library(grid)
grob1 <- grid::grobTree(grid::textGrob("TD\nBetter", x=-.15, y=.94, just = "centre",
gp=grid::gpar(col="red", fontsize=13, fontface="italic")))
grob2 <- grid::grobTree(grid::textGrob("RSTD\nBetter", x=-.15, y=.04, just = "centre",
gp=grid::gpar(col="red", fontsize=13, fontface="italic")))
#plot the best fitting model as a function of AI
plot2 <- ggplot(widedf, aes(x = TDRS_AsymmIdx, y= TDRS_TD_AIC)) +
theme_bw()+
geom_point(alpha = .5) +
geom_hline(yintercept=0, col = "red") +
labs(x = "Asymmetry Index", y = "AIC Difference between \nRSTD and TD models") +
annotation_custom(grob1) +
annotation_custom(grob2) +
theme(text=element_text(family="sans"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
gt2 <- ggplot_gtable(ggplot2::ggplot_build(plot2))
gt2$layout$clip[gt2$layout$name == "panel"] <- "off"
grid.draw(gt2)
#Asymmetry Index by age
ggplot(data=widedf, aes(x=Age, y=TDRS_AsymmIdx)) + geom_point() +
stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black") +
scale_y_continuous(limits = c(-1, 1)) +
theme_bw() +
ylab("Asymmetry Index") +
geom_hline(yintercept=0, linetype="dashed") +
theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank())
#linear age
lmAIAge <- lm(TDRS_AsymmIdx ~ Age_Z, data = widedf)
#quadratic age
lmAIAgeSq <- lm(TDRS_AsymmIdx ~ Age_Z+I(Age_Z^2), data = widedf)
anova(lmAIAge,lmAIAgeSq)
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 60 | 14.95345 | NA | NA | NA | NA |
| 59 | 13.59786 | 1 | 1.355589 | 5.881791 | 0.0183774 |
#quadratic better
lmReport(lmAIAgeSq)
| TDRS_AsymmIdx | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | -0.38 | -0.57 – -0.20 | -4.20 | 59.00 | <0.001 |
| Age_Z | -0.06 | -0.18 – 0.06 | -0.97 | 59.00 | 0.338 |
| Age_Z^2 | 0.17 | 0.03 – 0.31 | 2.43 | 59.00 | 0.018 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.108 / 0.078 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.02 | 0.95 | 0 | 0.16 |
| I(Age_Z^2) | 0.10 | 0.95 | 0 | 0.33 |
#AI by age - 2 lines approach
lmAIByAge <- twolines(TDRS_AsymmIdx ~ Age,data=widedf)
mean(widedf$TDRS_AsymmIdx)
[1] -0.2190999
sd(widedf$TDRS_AsymmIdx)
[1] 0.4998784
#Alpha+ by age
ggplot(data=widedf, aes(x=Age, y=TDRS_AlphaPos)) + geom_point() +
scale_y_continuous(limits = c(0, 1)) +
theme_bw() +
ylab(expression(paste(alpha,"+")))+
theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black")
#linear age
lmAPosAge <- lm(TDRS_AlphaPos ~ Age_Z, data = widedf)
#quadratic age
lmAPosAgeSq <- lm(TDRS_AlphaPos ~ Age_Z+I(Age_Z^2), data = widedf)
anova(lmAPosAge,lmAPosAgeSq)
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 60 | 2.120309 | NA | NA | NA | NA |
| 59 | 2.072704 | 1 | 0.0476044 | 1.355071 | 0.2490794 |
#quadratic not better
lmReport(lmAPosAge)
| TDRS_AlphaPos | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 0.21 | 0.16 – 0.26 | 8.72 | 60.00 | <0.001 |
| Age_Z | -0.02 | -0.07 – 0.03 | -0.85 | 60.00 | 0.401 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.012 / -0.005 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.01 | 0.95 | 0 | 0.13 |
#Alpha + by age - 2 lines approach
lmAPosByAge <- twolines(TDRS_AlphaPos ~ Age,data=widedf)
#Alpha- by age
ggplot(data=widedf, aes(x=Age, y=TDRS_AlphaNeg)) + geom_point() + scale_y_continuous(limits = c(0, 1)) +
theme_bw() +
ylab(expression(paste(alpha,"-")))+
theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + stat_smooth(method=lm,formula = y ~ x + I(x^2),color="black")
#linear age
lmAlphanegAge <- lm(TDRS_AlphaNeg ~ Age_Z, data = widedf)
#quadratic age
lmAlphanegAgeSq <- lm(TDRS_AlphaNeg ~ Age_Z+I(Age_Z^2), data = widedf)
anova(lmAlphanegAge,lmAlphanegAgeSq)
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 60 | 2.931939 | NA | NA | NA | NA |
| 59 | 2.542451 | 1 | 0.3894882 | 9.038443 | 0.0038793 |
lmReport(lmAlphanegAgeSq)
| TDRS_AlphaNeg | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 0.43 | 0.35 – 0.51 | 10.94 | 59.00 | <0.001 |
| Age_Z | -0.00 | -0.06 – 0.05 | -0.17 | 59.00 | 0.869 |
| Age_Z^2 | -0.09 | -0.15 – -0.03 | -3.01 | 59.00 | 0.004 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.133 / 0.103 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.00 | 0.95 | 0.00 | 0.00 |
| I(Age_Z^2) | 0.15 | 0.95 | 0.02 | 0.43 |
#Alpha - by age - 2 lines approach
lmAnegByAge <- twolines(TDRS_AlphaNeg ~ Age,data=widedf)
#Hit rate
mean(widedf$Hit)
[1] 0.5358259
sd(widedf$Hit)
[1] 0.1412183
#False alarm rate
mean(widedf$FA)
[1] 0.2357636
sd(widedf$FA)
[1] 0.1493429
#dprime
mean(widedf$dPrime)
[1] 0.9321024
sd(widedf$dPrime)
[1] 0.4827073
RiskySafeHits <- MemDF[,c("SubjectNumber","RespOld","AnyRisk")] %>%
group_by(SubjectNumber, AnyRisk) %>%
summarize(Hits = mean(RespOld))
## `summarise()` regrouping output by 'SubjectNumber' (override with `.groups` argument)
RiskySafeHits$AnyRiskNum <- RiskySafeHits$AnyRisk
RiskySafeHits$AnyRisk <- factor(RiskySafeHits$AnyRisk, levels = c(0,1), labels = c("Safe","Risky"))
RiskySafeHitswide <- tidyr::spread(RiskySafeHits[,c("SubjectNumber","Hits","AnyRisk")], AnyRisk, Hits)
#do the distributions deviate significantly from normality?
pander(shapiro.test(RiskySafeHitswide$Risky))
| Test statistic | P value |
|---|---|
| 0.9881 | 0.8097 |
pander(shapiro.test(RiskySafeHitswide$Safe))
| Test statistic | P value |
|---|---|
| 0.9807 | 0.438 |
#no, they don't (even with outliers)
pander(mean(RiskySafeHitswide$Risky))
0.5554
pander(sd(RiskySafeHitswide$Risky))
0.15
pander(mean(RiskySafeHitswide$Safe))
0.5225
pander(sd(RiskySafeHitswide$Safe))
0.1453
#paired t-test
RiskyVsSafeTOuts <- t.test(RiskySafeHitswide$Risky,RiskySafeHitswide$Safe,paired=TRUE)
pander(RiskyVsSafeTOuts)
| Test statistic | df | P value | Alternative hypothesis |
|---|---|---|---|
| 3.077 | 61 | 0.003126 * * | two.sided |
| mean of the differences |
|---|
| 0.03287 |
effectsize::cohens_d(RiskySafeHits$Hits,RiskySafeHits$AnyRisk, paired = TRUE)
| Cohens_d | CI | CI_low | CI_high |
|---|---|---|---|
| -0.3908301 | 0.95 | -0.6529867 | -0.1320529 |
#Hits
ggplot(widedf, aes(x = Age, y = Hit)) + stat_smooth(method=lm,formula = y ~ x ,se=TRUE,color="black") +geom_point() +ylab("Hits")
lmhitsage <- lm(Hit ~ Age_Z, data = widedf)
lmhitsagesq <- lm(Hit ~ Age_Z +I(Age_Z^2), data = widedf)
anova(lmhitsage, lmhitsagesq)
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 60 | 1.173963 | NA | NA | NA | NA |
| 59 | 1.170697 | 1 | 0.003266 | 0.1646 | 0.6864238 |
lmReport(lmhitsage)
| Hit | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 0.54 | 0.50 – 0.57 | 30.16 | 60.00 | <0.001 |
| Age_Z | 0.03 | -0.01 – 0.06 | 1.47 | 60.00 | 0.146 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.035 / 0.019 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.04 | 0.95 | 0 | 0.2 |
#False Alarms
ggplot(widedf, aes(x = Age, y = FA)) + stat_smooth(method=lm,formula = y ~ x ,se=TRUE,color="black") +geom_point() + ylab("False Alarms")
lmFAsage <- lm(FA ~ Age_Z, data = widedf)
lmFAsagesq <- lm(FA ~ Age_Z+I(Age_Z^2), data = widedf)
anova(lmFAsage, lmFAsagesq)
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 60 | 1.264663 | NA | NA | NA | NA |
| 59 | 1.240174 | 1 | 0.0244897 | 1.165074 | 0.2848087 |
lmReport(lmFAsage)
| FA | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 0.24 | 0.20 – 0.27 | 12.79 | 60.00 | <0.001 |
| Age_Z | 0.04 | 0.00 – 0.08 | 2.13 | 60.00 | 0.037 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.070 / 0.055 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.08 | 0.95 | 0 | 0.28 |
#D-Prime
ggplot(widedf, aes(x = Age, y = dPrime)) + stat_smooth(method=lm,formula = y ~ x ,se=TRUE,color="black") +geom_point() + ylab("d'")
lmdprimesage <- lm(dPrime ~ Age_Z, data = widedf)
lmReport(lmdprimesage)
| dPrime | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 0.93 | 0.81 – 1.05 | 15.50 | 60.00 | <0.001 |
| Age_Z | -0.11 | -0.23 – 0.01 | -1.84 | 60.00 | 0.070 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.054 / 0.038 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.06 | 0.95 | 0 | 0.24 |
lmdprimesagesq <- lm(dPrime ~ Age_Z+I(Age_Z^2), data = widedf)
lmReport(lmdprimesagesq)
| dPrime | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 0.93 | 0.75 – 1.11 | 10.27 | 59.00 | <0.001 |
| Age_Z | -0.11 | -0.23 – 0.01 | -1.83 | 59.00 | 0.073 |
| Age_Z^2 | -0.00 | -0.14 – 0.14 | -0.01 | 59.00 | 0.990 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.054 / 0.022 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.06 | 0.95 | 0 | 0.25 |
| I(Age_Z^2) | 0.00 | 0.95 | 0 | 0.00 |
anova(lmdprimesage,lmdprimesagesq)
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 60 | 13.45155 | NA | NA | NA | NA |
| 59 | 13.45151 | 1 | 3.85e-05 | 0.0001687 | 0.9896811 |
#make a new df including false alarms
MemDF2 <- merge(MemDF,widedf[,c("SubjectNumber","FA")],by.x="SubjectNumber", by.y = "SubjectNumber")
#also add confidence data; high-confidence responses are coded as 5 or 8; 6 or 7 are low confidence
MemDF2$Conf <- ifelse(MemDF2$MemResp==5|MemDF2$MemResp==8,1,0)
MemDF2$PositiveRPEC <- ifelse(MemDF2$PositiveRPE == "Positive Prediction Error", 1, ifelse(MemDF2$PositiveRPE == "Negative Prediction Error", -1, NA))
MemDF2$AbsRPEScale <- scale(MemDF2$AbsRPE)
MemDF2$FAScale <- scale(MemDF2$FA)
MemDF2$AIScale <- scale(MemDF2$AsymmIdx)
lmUnsignedRPE <- glmer(RespOld ~ AbsRPEScale+MemIdxScaled+Age_Z+I(Age_Z^2) +FAScale+ (1+AbsRPEScale+MemIdxScaled | SubjectNumber), family=binomial, data=MemDF2,control = glmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=1e6)))
glmerReport(lmUnsignedRPE)
| RespOld | ||||
|---|---|---|---|---|
| Predictors | z | p | Odds Ratios | CI |
| (Intercept) | 1.68 | 0.093 | 1.18 | 0.97 – 1.44 |
| AbsRPEScale | 3.04 | 0.002 | 1.08 | 1.03 – 1.13 |
| MemIdxScaled | -5.84 | <0.001 | 0.82 | 0.76 – 0.87 |
| Age_Z | 0.23 | 0.817 | 1.02 | 0.88 – 1.17 |
| Age_Z^2 | -0.23 | 0.821 | 0.98 | 0.85 – 1.14 |
| FAScale | 5.13 | <0.001 | 1.45 | 1.26 – 1.67 |
| Random Effects | ||||
| σ2 | 3.29 | |||
| τ00 SubjectNumber | 0.25 | |||
| τ11 SubjectNumber.AbsRPEScale | 0.01 | |||
| τ11 SubjectNumber.MemIdxScaled | 0.05 | |||
| ρ01 | 0.36 | |||
| -0.09 | ||||
| ICC | 0.09 | |||
| N SubjectNumber | 62 | |||
| Observations | 11309 | |||
| Marginal R2 / Conditional R2 | 0.048 / 0.131 | |||
lmUnsignedPosRPE <- glmer(RespOld ~ AbsRPEScale*PositiveRPEC+MemIdxScaled+Age_Z+I(Age_Z^2) +FAScale+ (1+AbsRPEScale*PositiveRPEC+MemIdxScaled || SubjectNumber), family=binomial, data=MemDF2,control = glmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=1e6)))
#This is the model that converges
lmAIRPEMaxConv <- glmer(RespOld ~ AIScale*AbsRPEScale*PositiveRPEC+MemIdxScaled+Age_Z+I(Age_Z^2) +FAScale+ (1+AbsRPEScale+MemIdxScaled || SubjectNumber), family=binomial, data=MemDF2,control = glmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=1e6)))
glmerReport(lmAIRPEMaxConv)
| RespOld | ||||
|---|---|---|---|---|
| Predictors | z | p | Odds Ratios | CI |
| (Intercept) | 1.45 | 0.147 | 1.17 | 0.94 – 1.46 |
| AIScale | -0.68 | 0.498 | 0.95 | 0.82 – 1.10 |
| AbsRPEScale | 4.75 | <0.001 | 1.19 | 1.11 – 1.28 |
| PositiveRPEC | -1.61 | 0.108 | 0.93 | 0.86 – 1.02 |
| MemIdxScaled | -5.83 | <0.001 | 0.82 | 0.76 – 0.87 |
| Age_Z | 0.32 | 0.750 | 1.02 | 0.89 – 1.17 |
| Age_Z^2 | -0.18 | 0.856 | 0.99 | 0.84 – 1.15 |
| FAScale | 4.86 | <0.001 | 1.41 | 1.23 – 1.61 |
| AIScale * AbsRPEScale | 0.43 | 0.664 | 1.01 | 0.96 – 1.07 |
| AIScale * PositiveRPEC | 1.94 | 0.052 | 1.07 | 1.00 – 1.15 |
|
AbsRPEScale * PositiveRPEC |
-0.15 | 0.878 | 0.99 | 0.93 – 1.07 |
|
(AIScale * AbsRPEScale) * PositiveRPEC |
3.45 | 0.001 | 1.12 | 1.05 – 1.19 |
| Random Effects | ||||
| σ2 | 3.29 | |||
| τ00 SubjectNumber | 0.25 | |||
| τ00 SubjectNumber.1 | 0.01 | |||
| τ00 SubjectNumber.2 | 0.05 | |||
| ICC | 0.07 | |||
| N SubjectNumber | 62 | |||
| Observations | 11309 | |||
| Marginal R2 / Conditional R2 | 0.049 / 0.115 | |||
#Is the model with AI significantly better than the one without?
anova(lmUnsignedRPE,lmAIRPEMaxConv)
| npar | AIC | BIC | logLik | deviance | Chisq | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|---|---|
| lmUnsignedRPE | 12 | 14766.53 | 14854.53 | -7371.264 | 14742.53 | NA | NA | NA |
| lmAIRPEMaxConv | 15 | 14752.02 | 14862.02 | -7361.011 | 14722.02 | 20.50578 | 3 | 0.0001333 |
#Plot
plotlabels <- c("Asymmetry Index (AI)", "PE Magnitude","PE Valence","Memory Trial Number","Linear Age","Quadratic Age","False Alarm Rate", "AI:PE Magnitude","AI:PE Valence","PE Magnitude:PE Valence","AI:PEMagnitude:PEValence")
names(plotlabels) <- c("AIScale", "AbsRPEScale","PositiveRPEC","MemIdxScaled","Age_Z","I(Age_Z^2)","FAScale","AIScale:AbsRPEScale","AIScale:PositiveRPEC","AbsRPEScale:PositiveRPEC","AIScale:AbsRPEScale:PositiveRPEC")
threewayintxplotPaper <- plot_model(lmAIRPEMaxConv, colors = "bw", show.values = TRUE,
value.offset = .4, order.terms=c(5,6,4,7,1,2,3,8,9,10,11),
title = "Fixed Effects",axis.labels = plotlabels,
vline.color = "grey", axis.lim = c(.3,3))
threewayintxplotPaper + ylab("Odds of Correct Memory Response") +
theme(text=element_text(family="sans",size=12),panel.background = element_rect(fill = "white", colour = "white"), panel.border = element_rect(colour = "black", fill=NA)) +
scale_y_log10(limits = c(.5,2))
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
#saved manually - 500x400 pixels
#The model that converges, but without scaled variables so I can plot the 3-way interaction (because this is the highest level interaction, we can use the unscaled numbers)
#also for some reason theme elements (grid color, background color, didn't work for this plot and the one above so i had to reset all of them here)
lmAIRPEMaxConvPlot <- glmer(RespOld ~ AsymmIdx*AbsRPE*PositiveRPE+MemIdxScaled+Age_Z+I(Age_Z^2) +FAScale+ (1+AbsRPE+MemIdxScaled || SubjectNumber), family=binomial, data=MemDF2,control = glmerControl(optimizer = "bobyqa",optCtrl=list(maxfun=1e6)))
AIRPEFig <- plot_model(lmAIRPEMaxConvPlot,
type = "pred",
terms = c("AbsRPE [all]", "AsymmIdx [-.8,0,.8]", "PositiveRPE"),
colors = c("#bdd7e7","#6baed6","#2171b5")) +
scale_y_continuous(breaks=c(.4,.5,.6,.7,.8,.9), labels = c("40%","50%","60%","70%","80%","90%"))
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
AIRPEFig <- AIRPEFig + theme_bw() + theme(
text=element_text(family="sans",size=12),
panel.background = element_rect(fill = "white", colour = "white"),
panel.border = element_rect(colour = "black", fill=NA),
panel.spacing = unit(1, "lines"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
AIRPEFig <- AIRPEFig + labs(title="Predicted Memory Accuracy by Choice Context",
x="Absolute Value of Prediction Error",
y= "Estimated Marginal Means for \n% Correct Memory Responses",
linetype="Asymmetry\nIndex",
color = "Asymmetry\nIndex")
AIRPEFig
#ggsave('plots/fig4b.png', width = 6, height = 3)
#saved manually - 600x400 pixels
#originally dospert was coded from 0 to 1. most papers report from 1 to 7, so rescaling accordingly
widedf$DOSPERT_rs <- widedf$DOSPERT*7
#linear age
lmDOSPERT <- lm(DOSPERT_rs ~ Age_Z, data = widedf)
#quadratic age
lmDOSPERTAgesq <- lm(DOSPERT_rs ~ Age_Z + I(Age_Z^2), data = widedf)
anova(lmDOSPERT,lmDOSPERTAgesq)
| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 60 | 60.65872 | NA | NA | NA | NA |
| 59 | 52.20750 | 1 | 8.451218 | 9.55077 | 0.0030481 |
#quadratic is better
lmReport(lmDOSPERTAgesq)
| DOSPERT_rs | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 4.17 | 3.81 – 4.53 | 23.31 | 59.00 | <0.001 |
| Age_Z | 0.15 | -0.09 – 0.39 | 1.27 | 59.00 | 0.208 |
| Age_Z^2 | -0.42 | -0.69 – -0.15 | -3.09 | 59.00 | 0.003 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.166 / 0.137 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.04 | 0.95 | 0.00 | 0.20 |
| I(Age_Z^2) | 0.16 | 0.95 | 0.02 | 0.44 |
DospertAge2l <- twolines(DOSPERT_rs ~ Age, data = widedf)
#plot
ggplot(data=widedf, aes(x=Age, y=DOSPERT_rs)) +
stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black") +
geom_point() +
theme_bw() +
ylab("Mean DOSPERT Score") +
theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank())
#risk predicting dospert
lmDOSPERTRisk <- lm(DOSPERT_rs ~ MeanRiskAll, data = widedf)
lmReport(lmDOSPERTRisk)
| DOSPERT_rs | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 4.01 | 3.42 – 4.60 | 13.58 | 60.00 | <0.001 |
| MeanRiskAll | -0.58 | -1.80 – 0.64 | -0.95 | 60.00 | 0.347 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.015 / -0.002 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| MeanRiskAll | 0.01 | 0.95 | 0 | 0.14 |
#correlation between dospert and task risk taking
pander(cor.test(widedf$DOSPERT_rs, widedf$MeanRiskAll))
| Test statistic | df | P value | Alternative hypothesis | cor |
|---|---|---|---|---|
| -0.9478 | 60 | 0.347 | two.sided | -0.1215 |
#95% CI of correlation
pander(cor.test(widedf$DOSPERT_rs, widedf$MeanRiskAll)$conf.int)
-0.3603 and 0.1323
#plot correlation
ggplot(data=widedf, aes(x=MeanRiskAll, y=DOSPERT_rs)) +
stat_smooth(method=lm,formula = y ~ x ,se=FALSE,color="grey",linetype=2) +
geom_point() +
theme_bw() +
xlab("Point Machine Risk Taking") +
ylab("Self-Reported Risk Taking") +
theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank())
ggplot(data=PPC, aes(x=Age, y=MeanRiskOverallEqEV_sim)) +
geom_hline(yintercept=0.5, linetype="dashed") +
stat_smooth(method=lm,formula = y ~ x + I(x^2),se=TRUE,color="black") +
geom_point() +
scale_y_continuous(limits = c(0, 1)) +
theme_bw() +
xlab("Real Participant's Age") +
ylab("Proportion Probabilistic Choices \n Equal-EV Trials - Simulated") +
theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank())
PPC_Reg <- lm(MeanRiskOverallEqEV_sim ~ Age_Z + I(Age_Z^2), data = PPC)
lmReport(PPC_Reg)
| MeanRiskOverallEqEV_sim | |||||
|---|---|---|---|---|---|
| Predictors | Estimates | CI | t | df | p |
| (Intercept) | 0.35 | 0.28 – 0.41 | 10.70 | 59.00 | <0.001 |
| Age_Z | -0.02 | -0.07 – 0.02 | -1.02 | 59.00 | 0.310 |
| Age_Z^2 | 0.05 | 0.00 – 0.10 | 2.16 | 59.00 | 0.035 |
| Observations | 62 | ||||
| R2 / R2 adjusted | 0.093 / 0.062 | ||||
| Parameter | Cohensf2Partial | f2_CI | f2_CI_low | f2_CI_high |
|---|---|---|---|---|
| Age_Z | 0.02 | 0.95 | 0 | 0.17 |
| I(Age_Z^2) | 0.08 | 0.95 | 0 | 0.29 |
PPC_RealData <- merge(PPC[,c("SubjectNumber","Age", "MeanRiskOverallEqEV_sim")], widedf[,c("SubjectNumber","MeanRiskOverallEqEV")])
#correlation between risk taking in simulated vs real data
pander(cor.test(PPC_RealData$MeanRiskOverallEqEV_sim,PPC_RealData$MeanRiskOverallEqEV))
| Test statistic | df | P value | Alternative hypothesis | cor |
|---|---|---|---|---|
| 0.2451 | 32 | 0.8079 | two.sided | 0.04329 |
#95% CI of correlation
pander(cor.test(PPC_RealData$MeanRiskOverallEqEV_sim,PPC_RealData$MeanRiskOverallEqEV)$conf.int)
-0.2993 and 0.376
ggplot(data=PPC_RealData, aes(x=MeanRiskOverallEqEV, y=MeanRiskOverallEqEV_sim)) +
stat_smooth(method=lm,formula = y ~ x ,se=TRUE,color="black") +
geom_point() +
scale_x_continuous(limits = c(0, 1)) +
scale_y_continuous(limits = c(0, 1)) +
theme_bw() +
xlab("Proportion Equal-EV \nProbabilistic Choices - Real") +
ylab("Proportion Equal-EV \nProbabilistic Choices - Simulated") +
theme(text=element_text(family="sans",size=16),panel.grid.major = element_blank(), panel.grid.minor = element_blank())
#compute AI for simulated data
TDRSRecovDF$AsymmIdx <- (TDRSRecovDF$RealAlphaPos-TDRSRecovDF$RealAlphaNeg)/(TDRSRecovDF$RealAlphaPos+TDRSRecovDF$RealAlphaNeg)
# make variables comparing TD to TDRS model fit. higher numbers mean better fit for TD model, while lower numbers mean better fit for TDRS model
#compare AIC for data generated by TDRS
TDRecovDF$TDRSvsTD_AIC <- -TDRecovDF$AIC_TD+TDRecovDF$AIC_TDRS
#compare AIC for data generated by TDRS
TDRSRecovDF$TDRSvsTD_AIC <- -TDRSRecovDF$AIC_TD+TDRSRecovDF$AIC_TDRS
#TD model - plot the proportion of sims best fit by the TD model vs. TDRS
PropBestFitTD <- nrow(subset(TDRecovDF,TDRSvsTD_AIC>0))/nrow(TDRecovDF)
str <- paste("proportion \nof TD subs\n best fit by \nTD model:", toString(PropBestFitTD))
AICTDDiffPlot <- ggplot(TDRecovDF, aes(TDRSvsTD_AIC)) + geom_histogram() + theme_bw() + labs(x = "AIC difference (Higher = TD is Better)") + geom_vline(xintercept=0) + annotate("text", x=-100, y=5000,label=str) +theme_bw()
AICTDDiffPlot
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#TDRS model - plot the proportion of sims best fit by the TDRS model vs. TD
PropBestFitTDRS <- nrow(subset(TDRSRecovDF,TDRSvsTD_AIC<0))/nrow(TDRSRecovDF)
str <- paste("proportion \nof TDRS subs\n best fit by \nTDRS model:", toString(PropBestFitTDRS))
AICTDRSDiffPlot <- ggplot(TDRSRecovDF, aes(TDRSvsTD_AIC)) + geom_histogram() + theme_bw() + labs(x = "AIC difference (Lower = TDRS is Better)") + geom_vline(xintercept=0) + annotate("text", x=120, y=1700,label=str) +theme_bw()
AICTDRSDiffPlot
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
grob1 <- grobTree(textGrob("TD\nBetter", x=-.15, y=.94, just = "centre",
gp=gpar(col="red", fontsize=13, fontface="italic")))
grob2 <- grobTree(textGrob("RSTD\nBetter", x=-.15, y=.04, just = "centre",
gp=gpar(col="red", fontsize=13, fontface="italic")))
#plot the best fitting model as a function of AI
plot <- ggplot(TDRSRecovDF, aes(x = AsymmIdx, y= TDRSvsTD_AIC)) +
geom_point(alpha = .2) +
geom_hline(yintercept=0, col = "red") +
labs(x = "Asymmetry Index", y = "AIC Difference between \nRSTD and TD models") +
annotation_custom(grob1) +
annotation_custom(grob2)
gt <- ggplot_gtable(ggplot_build(plot))
gt$layout$clip[gt$layout$name == "panel"] <- "off"
grid.draw(gt)
#ggsave('plots/figS1.png', width = 4, height = 3)
#simulated subs with a more extreme AIs are better fit by TD, while those with less extreme AIs show less of a difference in model fitting (which makes sense)
#real vs. recovered alpha
acor <- toString(round(cor(TDRecovDF$RealAlpha,TDRecovDF$RecAlpha),2))
acorann <- paste('r =',acor)
ggplot(TDRecovDF, aes(RealAlpha,RecAlpha)) + geom_point() + geom_abline(intercept = 0, slope = 1,color="blue") +annotate("text", x=.25, y=.9,label=acorann) +theme_bw()
#real vs. recovered beta
bcor <- toString(round(cor(TDRecovDF$RealBeta,TDRecovDF$RecBeta),2))
bcorann = paste('r =',bcor)
ggplot(TDRecovDF, aes(RealBeta,RecBeta)) + geom_point()+ geom_abline(intercept = 0, slope = 1,color="blue") +annotate("text", x=5, y=20,label=bcorann) +theme_bw()
#bin recovery for plotting later
TDRSRecovDF$alphaposlevel <- ceiling(TDRSRecovDF$RealAlphaPos/.05)*.05
TDRSRecovDF$alphaneglevel <- ceiling(TDRSRecovDF$RealAlphaNeg/.05)*.05
TDRSRecovDF$betalevel <- ceiling(TDRSRecovDF$RealBeta)
#correlation between real and recovered alpha positive
aposcor <- toString(round(cor(TDRSRecovDF$RealAlphaPos,TDRSRecovDF$RecAlphaPos),2))
aposcorann = paste('r =',aposcor)
ggplot(TDRSRecovDF, aes(RealAlphaPos,RecAlphaPos)) + geom_point() + annotate("text", x=.25, y=.9,label=aposcorann)+ geom_abline(intercept = 0, slope = 1,color="blue") +theme_bw()
#correlation between real and recovered alpha negative
anegcor <- toString(round(cor(TDRSRecovDF$RealAlphaNeg,TDRSRecovDF$RecAlphaNeg),2))
anegcorann = paste('r =',anegcor)
ggplot(TDRSRecovDF, aes(RealAlphaNeg,RecAlphaNeg)) + geom_point()+ annotate("text", x=.25, y=.9,label=anegcorann)+ geom_abline(intercept = 0, slope = 1,color="blue") +theme_bw()
#correlation between real and recovered beta
bcorTDRS <- toString(round(cor(TDRSRecovDF$RealBeta,TDRSRecovDF$RecBeta),2))
bcorTDRSann <- paste('r =',bcorTDRS)
ggplot(TDRSRecovDF, aes(RealBeta,RecBeta)) + geom_point() + geom_point()+annotate("text", x=5, y=20,label=bcorTDRSann) + geom_abline(intercept = 0, slope = 1,color="blue") +theme_bw()